#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "IComplexMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class ComplexHermitianMatrix : public IComplexMatrix {
	internal:
		doublecomplex * ap;
	private:
		long * ipiv;
		int m_rownum;
		int m_colnum;
		bool m_isFactorized;

		ComplexHermitianMatrix(ComplexHermitianMatrix ^ orig) {
			m_rownum = orig->m_rownum;
			m_colnum = orig->m_colnum;

			int size = m_rownum * (m_rownum + 1) / 2;
			ap = new doublecomplex[size];

			for (int i = 0; i < size; i++) {
				ap[i] = orig->ap[i];
			}
			ipiv = new long[m_rownum];

			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}

			m_isFactorized = orig->m_isFactorized;

		};
	public:
		ComplexHermitianMatrix(int rownum, int colnum) {
			if (rownum != colnum) {
				throw gcnew LatoolException("Number of rows and columns must be equal for symmetric matrix.");
			}
			m_rownum = rownum;
			m_colnum = colnum;

			int size = rownum * (rownum + 1) / 2;
			ap = new doublecomplex[size];

			for (int i = 0; i < size; i++) {
				ap[i].i = 0.0;
				ap[i].r = 0.0;
			}
      ipiv = new long[rownum];

			m_isFactorized = false;
		};

		~ComplexHermitianMatrix(){
			delete[] ap;
		};

		virtual IComplexMatrix ^ Clone() {		
			return (IComplexMatrix^) gcnew ComplexHermitianMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property Complex default[int, int] {
			Complex get(int i, int j) {
				bool conj = false; 
				if (i > j) {
					int temp = i;
					i = j;
					j = temp;
					conj = true;
				}
				Complex ret = Complex(ap[i + j * (j + 1) / 2].r, ap[i + j * (j + 1) / 2].i);
				if (conj) {
					return Complex::Conj(ret);
				} else {
					return ret;
				}
			}
			void set(int i, int j, Complex value) {
				if (i > j) {
					int temp = i;
					i = j;
					j = temp;
					value = Complex::Conj(value);
				}
				doublecomplex c;
				c.r = value.Real;
				c.i = value.Imag;
				ap[i + j * (j + 1) / 2] = c;
				m_isFactorized = false;
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::ComplexHermitian; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		}

		virtual void Invert() {

			char uplo = 'U';
			long n = m_rownum;
			long info;

				
			doublecomplex * work = new doublecomplex[Math::Max(1, n)];

			zhptri_(&uplo, &n, ap, ipiv, work, &info);

			delete[] work;

			m_isFactorized = false;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("the i-th diagonal element of the Cholesky factor (and therefore the factor itself) is zero, and the inversion could not be completed.");
			}

		};

		virtual void Factorize() {

			char uplo = 'U';
			long n = m_rownum;
			long info;

			zhptrf_(&uplo, &n, ap, ipiv, &info);

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("The factorization has been completed, but D is exactly singular. Division by 0 will occur if you use D for solving a system of linear equations.");
			}

			m_isFactorized = true;

		};

		virtual void Solve(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			char uplo = 'U';
			long n = m_rownum;
			long nrhs = b->ColNum;
			long ldb = Math::Max(1, n);
			long info;

			zhptrs_(&uplo, &n, &nrhs, ap, ipiv, gb->dat, &ldb, &info);

			m_isFactorized = false;
			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			char uplo = 'U';
			long n = m_rownum;
			long nrhs = b->ColNum;
			long ldb = Math::Max(1, n);
			long info;

			zppsv_(&uplo, &n, &nrhs, ap, gb->dat, &ldb, &info);
			
			if (info == 0) {
				m_isFactorized = true;
				return;
			}

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			zhpsv_(&uplo, &n, &nrhs, ap, ipiv, gb->dat, &ldb, &info);

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("The factorization has been completed, but D is exactly singular, so the solution could not be computed.");
			}

			m_isFactorized = true;

		};

		void SolveHermitianEigenproblem(array<double>^ values, IMatrix^ vectors) {

			char jobz = 'V';
			char uplo = 'U';
			long n = m_rownum;
			doublecomplex * work = new doublecomplex[2 * n - 1];
			long ldz = n;
			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;
			doublecomplex * z = ((ComplexGeneralMatrix ^) vectors)->dat;

			double * rwork = new double[3 * n - 2];
			long info = 0;
		
			zhpev_(&jobz, &uplo, &n, ap, w, z, &ldz, work, rwork, &info);

			delete[] work;
			delete[] rwork;
			
			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException(String::Format("{i} elements of an intermediate tridiagonal form did not converge to zero.", info));
			}

		};

		void SolveHermitianEigenproblem(array<double>^ values) {

			char jobz = 'N';
			char uplo = 'U';
			long n = m_rownum;
			doublecomplex * work = new doublecomplex[2 * n - 1];
			long ldz = n;
			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;
			doublecomplex z;

			double * rwork = new double[3 * n - 2];
			long info = 0;
		
			zhpev_(&jobz, &uplo, &n, ap, w, &z, &ldz, work, rwork, &info);

			delete[] work;
			delete[] rwork;
			
			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException(String::Format("{i} elements of an intermediate tridiagonal form did not converge to zero.", info));
			}

		};

		void SolveGeneralizedHermitianDefiniteEigenproblem(IMatrix^ b, array<double>^ values, IMatrix^ vectors) {
			
			long itype = 1;
			char jobz = 'V';
			char uplo = 'U';
			long n = m_rownum;
			doublecomplex * bp = ((ComplexHermitianMatrix ^) b)->ap;

			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;

			doublecomplex * z = ((ComplexGeneralMatrix ^) vectors)->dat;
			long ldz = n;
			doublecomplex * work = new doublecomplex[2 * n - 1];
			double * rwork = new double[3 * n - 2];
			long info = 0;

			zhpgv_(&itype, &jobz, &uplo, &n, ap, bp, w, z, &ldz, work, rwork, &info);

			delete[] work;
			delete[] rwork;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("off-diagonal elements of an intermediate tridiagonal did not converge to zero");
			}

		};

		void SolveGeneralizedHermitianDefiniteEigenproblem(IMatrix^ b, array<double>^ values) {
			
			long itype = 1;
			char jobz = 'N';
			char uplo = 'U';
			long n = m_rownum;
			doublecomplex * bp = ((ComplexHermitianMatrix ^) b)->ap;

			pin_ptr<double> ptr =&(values[0]);
			double * w = ptr;

			doublecomplex z;
			long ldz = n;
			doublecomplex * work = new doublecomplex[2 * n - 1];
			double * rwork = new double[3 * n - 2];
			long info = 0;

			zhpgv_(&itype, &jobz, &uplo, &n, ap, bp, w, &z, &ldz, work, rwork, &info);

			delete[] work;
			delete[] rwork;

			if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			}

			if (info > 0) {
				throw gcnew LatoolException("off-diagonal elements of an intermediate tridiagonal did not converge to zero");
			}

		};

	};
}